Autocorrelación espacial

Autor

Omar E. Barrantes Sotela

Fecha de Publicación

14 de septiembre de 2022

Introducción

La autocorrelación espacial es un concepto importante en las estadísticas espaciales. Su fundamento teórico es lo que permite la interpolación espacial. En muchas ocasiones su calculo y propiedades son malentendidas.

La autocorrelación (ya sea espacial o no) es una medida de similitud (correlación) entre observaciones cercanas. Para comprender la autocorrelación espacial, primero es útil considerar la autocorrelación temporal.

Autocorrelación temporal

Si mide algo sobre el mismo objeto a lo largo del tiempo, por ejemplo, el peso de una persona, es probable que dos observaciones que estén cercanas entre sí en el tiempo también sean similares en la medición. Digamos que en un par de años su peso pasó de \(50\) a \(80\) kg. Es poco probable que fuera \(60\) kg un día, \(50\) kg el siguiente y \(80\) al día siguiente. Por el contrario, probablemente subió gradualmente, con la reducción gradual ocasional, o incluso revertir en la dirección. Lo mismo puede ser cierto con su cuenta bancaria, pero eso también puede tener una tendencia mensual marcada. Para medir el grado de asociación a lo largo del tiempo, podemos calcular la correlación de cada observación con la siguiente observación.

En este caso d es un vector de observaciones diarias.

set.seed(123)
d <- sample(100, 10)
d
 [1] 31 79 51 14 67 42 50 43 97 25

Cálculo de la auto-correlación

a <- d[-length(d)]
b <- d[-1]
plot(a, b, xlab='t', ylab='t-1')

cor(a, b) #Calculo de la correlación
[1] -0.55068

La autocorrelación calculada anteriormente es muy pequeña. Aunque esta es una muestra aleatoria, (casi) nunca obtiene un valor de cero. Calculamos la autocorrelación de “un retraso”, es decir, comparamos cada valor con su vecino inmediato y no con otros valores cercanos.

Después de ordenar los números en la variable d, la autocorrelación se vuelve muy fuerte (como era de esperar).

d <- sort(d) #ordena los valores
d
 [1] 14 25 31 42 43 50 51 67 79 97
a <- d[-length(d)]
b <- d[-1]
plot(a, b, xlab='t', ylab='t-1')

cor(a,b)
[1] 0.9711034

La función acf muestra la autocorrelación calculada de una manera ligeramente diferente para varios retrasos (es 1 para cada punto en sí mismo, muy alta cuando se compara con el vecino más cercano, y se reduce gradualmente).

acf(d)

Autocorrelación espacial

El concepto de autocorrelación espacial es una extensión de la autocorrelación temporal. Aunque es un poco más complicado. El tiempo es unidimensional, y solo va en una dirección, siempre adelante. Los objetos espaciales tienen (al menos) dos dimensiones y formas complejas, y puede no ser obvio cómo determinar qué está “cerca”.

Las medidas de autocorrelación espacial describen el grado en que las observaciones (valores) en ubicaciones espaciales (ya sean puntos, áreas o celdas de trama) son similares entre sí. Entonces necesitamos dos cosas: observaciones y ubicaciones.

La autocorrelación espacial en una variable puede ser exógena (es causada por otra variable espacialmente autocorrelacionada, por ejemplo, lluvia) o endógena (es causada por el proceso en juego, por ejemplo, la diseminación de una enfermedad).

Una estadística comúnmente utilizada que describe la autocorrelación espacial es la I de Moran, y discutiremos esto aquí en detalle. Otros índices incluyen C de Geary y, para datos binarios, el índice de recuento de uniones. El semi-variograma también expresa la cantidad de autocorrelación espacial en un conjunto de datos (el cual se estudiará en el tema de interpolación) (Barrantes-Sotela & Solano-Mayorga, 2020).

Demostración

library(rgdal)
library(raster)

#Carga un shapefile en R
shape <- "./shp/San_Pedro_CU.shp"
p <- readOGR(dsn = shape,layer = 'San_Pedro_CU')
OGR data source with driver: ESRI Shapefile 
Source: "G:\UNA\Clases\ECG-UNA\MC\gei411-mc\demo\demo03\shp\San_Pedro_CU.shp", layer: "San_Pedro_CU"
with 14 features
It has 28 fields
Integer64 fields read as strings:  ugm tot_viv tot_hog hog_nbi tot_pp pp_noCCSS 
print(proj4string(p))  # Proyeccion del shapefile
[1] "+proj=tmerc +lat_0=0 +lon_0=-84 +k=0.9999 +x_0=500000 +y_0=0 +ellps=WGS84 +units=m +no_defs"

Sí se supone el interés en la autocorrelación espacial en la variable “Necesidades Básicas Insatisfechas”. Si hubiera autocorrelación espacial, las Unidades Mínimas Geocensalescon valor similar estarían agrupadas espacialmente.

Se construye un diagrama de los polígonos. Se utiliza la función coordinates para obtener los centroides de los polígonos para colocar las etiquetas.

par(mai=c(0,0,0,0))
plot(p, col=2:7)
xy <- coordinates(p)
points(xy, cex=6, pch=20, col='white')
text(p, 'cod', cex=1.5)

Adyacencia de los polígonos

Ahora necesitamos determinar qué polígonos están “cerca” y cómo cuantificar eso. Aquí usaremos la adyacencia como criterio. Para encontrar polígonos adyacentes, podemos usar el paquete spdep.

library(knitr)
library(spdep)
w <- poly2nb(p, row.names=p$cod)
class(w)
[1] "nb"
summary(w)
Neighbour list object:
Number of regions: 14 
Number of nonzero links: 64 
Percentage nonzero weights: 32.65306 
Average number of links: 4.571429 
Link number distribution:

3 4 5 6 8 
4 4 2 3 1 
4 least connected regions:
O M A B with 3 links
1 most connected region:
K with 8 links

summary(w) genera un resumen de los vecinos. El número promedio de vecinos (polígonos adyacentes) es 4.57, 4 con 3 enlaces, 4 con 4 enlaces, 2 con 5 enlaces, 3 con 6 enlaces y 1 con 8 (¿Cúal es ese? Respuesta: k.).

Para mayores detalles se puede observar la estructura (str) de w.

str(w)
List of 14
 $ : int [1:6] 2 3 7 8 9 14
 $ : int [1:4] 1 7 8 14
 $ : int [1:4] 1 8 9 14
 $ : int [1:3] 5 7 8
 $ : int [1:5] 4 6 7 8 9
 $ : int [1:3] 5 8 9
 $ : int [1:5] 1 2 4 5 8
 $ : int [1:8] 1 2 3 4 5 6 7 9
 $ : int [1:6] 1 3 5 6 8 14
 $ : int [1:4] 11 12 13 14
 $ : int [1:4] 10 12 13 14
 $ : int [1:3] 10 11 13
 $ : int [1:3] 10 11 12
 $ : int [1:6] 1 2 3 9 10 11
 - attr(*, "class")= chr "nb"
 - attr(*, "region.id")= chr [1:14] "H" "I" "G" "O" ...
 - attr(*, "call")= language poly2nb(pl = p, row.names = p$cod)
 - attr(*, "type")= chr "queen"
 - attr(*, "sym")= logi TRUE

Gráfico de los enlaces entre polígonos

plot(p, col='gray', border='blue', lwd=2) +
plot(w, xy, col='red', lwd=2, add=TRUE)
integer(0)

Podemos transformar w en una matriz de ponderaciones espaciales. Una matriz de ponderaciones espaciales refleja la intensidad de la relación geográfica entre las observaciones.

wm <- nb2mat(w, style='B')
kable(wm)
Tabla 1: Matriz de ponderaciones espaciales
H 0 1 1 0 0 0 1 1 1 0 0 0 0 1
I 1 0 0 0 0 0 1 1 0 0 0 0 0 1
G 1 0 0 0 0 0 0 1 1 0 0 0 0 1
O 0 0 0 0 1 0 1 1 0 0 0 0 0 0
N 0 0 0 1 0 1 1 1 1 0 0 0 0 0
M 0 0 0 0 1 0 0 1 1 0 0 0 0 0
L 1 1 0 1 1 0 0 1 0 0 0 0 0 0
K 1 1 1 1 1 1 1 0 1 0 0 0 0 0
J 1 0 1 0 1 1 0 1 0 0 0 0 0 1
D 0 0 0 0 0 0 0 0 0 0 1 1 1 1
C 0 0 0 0 0 0 0 0 0 1 0 1 1 1
A 0 0 0 0 0 0 0 0 0 1 1 0 1 0
B 0 0 0 0 0 0 0 0 0 1 1 1 0 0
F 1 1 1 0 0 0 0 0 1 1 1 0 0 0

I’ Moran

El cálculo de la I’ Moran corresponde a la siguiente formula

\[ I=\frac{n} {\sum_{i=1}^{n} {(y_i - \bar{y})^2}} \frac{\sum_{i}\sum _{j}w_{ij}(y_{i}-{\bar {y}})(y_{j}-{\bar {y}})} {\sum_{i}\sum_{j}w_{ij}} \]

Otra forma de verlo sería

\[ I={\frac {n}{w}}{\frac {\sum _{i}\sum _{j}w_{ij}(y_{i}-{\bar {y}})(y_{j}-{\bar {y}})}{\sum _{i}(y_{i}-{\bar {y}})^{2}}} \]

I’ Moran es una versión expandida de la fórmula para calcular el coeficiente de correlación. Lo principal que se agregó es la matriz de ponderaciones espaciales.

Primero se determina el número de observaciones

n <- length(p)
n
[1] 14

Se obtiene los valores \(y\) y \(\bar{y}\)

y <- p$hog_NBI_p
ybar <- mean(y)

Ahora se obtiene \((y_i - \bar{y})(y_j - \bar{y})\), para todos los pares de datos. Se puede realizar con dos métodos.

Método 1

dy <- y - ybar
g <- expand.grid(dy, dy)
yiyj <- g[,1] * g[,2]

Método 2

yi <- rep(dy, each=n)
yj <- rep(dy)
yiyj <- yi * yj

Se procede a obtener una matriz de los pares multiplicados

pm <- matrix(yiyj, ncol=n)
pm
            [,1]        [,2]       [,3]        [,4]        [,5]        [,6]
 [1,]  24.566184  -24.799844  30.018256  -39.322180  -62.914780   72.247027
 [2,] -24.799844   25.035727 -30.303773   39.696191   63.513191  -72.934202
 [3,]  30.018256  -30.303773  36.680327  -48.049109  -76.877709   88.281098
 [4,] -39.322180   39.696191 -48.049109   62.941556  100.705356 -115.643137
 [5,] -62.914780   63.513191 -76.877709  100.705356  161.126756 -185.026937
 [6,]  72.247027  -72.934202  88.281098 -115.643137 -185.026937  212.472270
 [7,] -43.089066   43.498906 -52.651994   68.971070  110.352470 -126.721223
 [8,] -62.914780   63.513191 -76.877709  100.705356  161.126756 -185.026937
 [9,] -62.914780   63.513191 -76.877709  100.705356  161.126756 -185.026937
[10,] 102.282984 -103.255844 124.983056 -163.720580 -261.949980  300.805427
[11,]  43.301484  -43.713344  52.911556  -69.311080 -110.896480  127.345927
[12,]  24.566184  -24.799844  30.018256  -39.322180  -62.914780   72.247027
[13,]   1.023148   -1.032880   1.250220   -1.637716   -2.620316    3.008991
[14,]  -2.049837    2.069334  -2.504766    3.281098    5.249698   -6.028394
             [,7]        [,8]        [,9]       [,10]       [,11]      [,12]
 [1,]  -43.089066  -62.914780  -62.914780  102.282984   43.301484  24.566184
 [2,]   43.498906   63.513191   63.513191 -103.255844  -43.713344 -24.799844
 [3,]  -52.651994  -76.877709  -76.877709  124.983056   52.911556  30.018256
 [4,]   68.971070  100.705356  100.705356 -163.720580  -69.311080 -39.322180
 [5,]  110.352470  161.126756  161.126756 -261.949980 -110.896480 -62.914780
 [6,] -126.721223 -185.026937 -185.026937  300.805427  127.345927  72.247027
 [7,]   75.578184  110.352470  110.352470 -179.404266  -75.950766 -43.089066
 [8,]  110.352470  161.126756  161.126756 -261.949980 -110.896480 -62.914780
 [9,]  110.352470  161.126756  161.126756 -261.949980 -110.896480 -62.914780
[10,] -179.404266 -261.949980 -261.949980  425.862184  180.288684 102.282984
[11,]  -75.950766 -110.896480 -110.896480  180.288684   76.325184  43.301484
[12,]  -43.089066  -62.914780  -62.914780  102.282984   43.301484  24.566184
[13,]   -1.794602   -2.620316   -2.620316    4.259948    1.803448   1.023148
[14,]    3.595413    5.249698    5.249698   -8.534637   -3.613137  -2.049837
            [,13]       [,14]
 [1,]  1.02314847 -2.04983724
 [2,] -1.03288010  2.06933418
 [3,]  1.25021990 -2.50476582
 [4,] -1.63771582  3.28109847
 [5,] -2.62031582  5.24969847
 [6,]  3.00899133 -6.02839439
 [7,] -1.79460153  3.59541276
 [8,] -2.62031582  5.24969847
 [9,] -2.62031582  5.24969847
[10,]  4.25994847 -8.53463724
[11,]  1.80344847 -3.61313724
[12,]  1.02314847 -2.04983724
[13,]  0.04261276 -0.08537296
[14,] -0.08537296  0.17104133

Ahora es posible multiplicar esta matriz con los pesos para establecer en cero el valor para los pares que no son adyacentes.

pmw <- pm * wm
pmw
        [,1]       [,2]       [,3]      [,4]      [,5]      [,6]      [,7]
H   0.000000 -24.799844  30.018256   0.00000    0.0000    0.0000 -43.08907
I -24.799844   0.000000   0.000000   0.00000    0.0000    0.0000  43.49891
G  30.018256   0.000000   0.000000   0.00000    0.0000    0.0000   0.00000
O   0.000000   0.000000   0.000000   0.00000  100.7054    0.0000  68.97107
N   0.000000   0.000000   0.000000 100.70536    0.0000 -185.0269 110.35247
M   0.000000   0.000000   0.000000   0.00000 -185.0269    0.0000   0.00000
L -43.089066  43.498906   0.000000  68.97107  110.3525    0.0000   0.00000
K -62.914780  63.513191 -76.877709 100.70536  161.1268 -185.0269 110.35247
J -62.914780   0.000000 -76.877709   0.00000  161.1268 -185.0269   0.00000
D   0.000000   0.000000   0.000000   0.00000    0.0000    0.0000   0.00000
C   0.000000   0.000000   0.000000   0.00000    0.0000    0.0000   0.00000
A   0.000000   0.000000   0.000000   0.00000    0.0000    0.0000   0.00000
B   0.000000   0.000000   0.000000   0.00000    0.0000    0.0000   0.00000
F  -2.049837   2.069334  -2.504766   0.00000    0.0000    0.0000   0.00000
        [,8]        [,9]      [,10]      [,11]      [,12]    [,13]     [,14]
H  -62.91478  -62.914780   0.000000   0.000000   0.000000 0.000000 -2.049837
I   63.51319    0.000000   0.000000   0.000000   0.000000 0.000000  2.069334
G  -76.87771  -76.877709   0.000000   0.000000   0.000000 0.000000 -2.504766
O  100.70536    0.000000   0.000000   0.000000   0.000000 0.000000  0.000000
N  161.12676  161.126756   0.000000   0.000000   0.000000 0.000000  0.000000
M -185.02694 -185.026937   0.000000   0.000000   0.000000 0.000000  0.000000
L  110.35247    0.000000   0.000000   0.000000   0.000000 0.000000  0.000000
K    0.00000  161.126756   0.000000   0.000000   0.000000 0.000000  0.000000
J  161.12676    0.000000   0.000000   0.000000   0.000000 0.000000  5.249698
D    0.00000    0.000000   0.000000 180.288684 102.282984 4.259948 -8.534637
C    0.00000    0.000000 180.288684   0.000000  43.301484 1.803448 -3.613137
A    0.00000    0.000000 102.282984  43.301484   0.000000 1.023148  0.000000
B    0.00000    0.000000   4.259948   1.803448   1.023148 0.000000  0.000000
F    0.00000    5.249698  -8.534637  -3.613137   0.000000 0.000000  0.000000
attr(,"call")
nb2mat(neighbours = w, style = "B")

Ahora se suman los valores, para obtener un aproximado del I de Moran:

\[ \sum _{i}\sum _{j}w_{ij}(y_{i}-{\bar {y}})(y_{j}-{\bar {y}}) \]

spmw <- sum(pmw)
spmw
[1] 1065.038

El siguiente paso es dividir este valor por la suma de los pesos.

smw <- sum(wm)
sw  <- spmw / smw

Se calcula la varianza inversa de y

vr <- n / sum(dy^2)

El último paso para calcular el I de Moran.

MI <- vr * sw
MI
[1] 0.1609378

Valor esperado de I’ Moran

Esta es una forma simple (pero cruda) de estimar el valor esperado de I. de Moran, es decir, el valor que obtendría en ausencia de autocorrelación espacial (si los datos fueran espacialmente aleatorios). Por supuesto, nunca se espera eso, pero así es como se hace en las estadísticas. Se debe tener en cuenta que el valor esperado se acerca a cero si n se vuelve grande, pero que no lo suficiente para valores pequeños de n.

EI <- -1/(n-1)
EI
[1] -0.07692308

Librerias para el cálculo del I’ Moran

Después de hacer esto a mano, ahora se usará el paquete spdep para calcular el I de Moran y realizar una prueba de significancia. Para hacer esto, se necesita crear un objeto de pesos espaciales tipo listw (en lugar de la matriz que utilizamos anteriormente). Para obtener el mismo valor que el anterior, usamos style = 'B' para usar pesos de distancia binarios (TRUE / FALSE).

ww <-  nb2listw(w, style='B')
ww
Characteristics of weights list object:
Neighbour list object:
Number of regions: 14 
Number of nonzero links: 64 
Percentage nonzero weights: 32.65306 
Average number of links: 4.571429 

Weights style: B 
Weights constants summary:
   n  nn S0  S1   S2
B 14 196 64 128 1288

Ahora podemos usar la función moran. Echa un vistazo con ?moran. La función se define como moran (y, ww, n, Szero (ww)). Tenga en cuenta los argumentos impares n y S0. Creo que son raros, porque “ww” tiene esa información. De todos modos, los suministramos y funciona. Probablemente haya casos en los que tenga sentido usar otros valores.

moran(p$hog_NBI_p, ww, n=length(ww$neighbours), S0=Szero(ww))
$I
[1] 0.1609378

$K
[1] 2.158294
Szero(ww)
[1] 64
# is the same as
pmw
        [,1]       [,2]       [,3]      [,4]      [,5]      [,6]      [,7]
H   0.000000 -24.799844  30.018256   0.00000    0.0000    0.0000 -43.08907
I -24.799844   0.000000   0.000000   0.00000    0.0000    0.0000  43.49891
G  30.018256   0.000000   0.000000   0.00000    0.0000    0.0000   0.00000
O   0.000000   0.000000   0.000000   0.00000  100.7054    0.0000  68.97107
N   0.000000   0.000000   0.000000 100.70536    0.0000 -185.0269 110.35247
M   0.000000   0.000000   0.000000   0.00000 -185.0269    0.0000   0.00000
L -43.089066  43.498906   0.000000  68.97107  110.3525    0.0000   0.00000
K -62.914780  63.513191 -76.877709 100.70536  161.1268 -185.0269 110.35247
J -62.914780   0.000000 -76.877709   0.00000  161.1268 -185.0269   0.00000
D   0.000000   0.000000   0.000000   0.00000    0.0000    0.0000   0.00000
C   0.000000   0.000000   0.000000   0.00000    0.0000    0.0000   0.00000
A   0.000000   0.000000   0.000000   0.00000    0.0000    0.0000   0.00000
B   0.000000   0.000000   0.000000   0.00000    0.0000    0.0000   0.00000
F  -2.049837   2.069334  -2.504766   0.00000    0.0000    0.0000   0.00000
        [,8]        [,9]      [,10]      [,11]      [,12]    [,13]     [,14]
H  -62.91478  -62.914780   0.000000   0.000000   0.000000 0.000000 -2.049837
I   63.51319    0.000000   0.000000   0.000000   0.000000 0.000000  2.069334
G  -76.87771  -76.877709   0.000000   0.000000   0.000000 0.000000 -2.504766
O  100.70536    0.000000   0.000000   0.000000   0.000000 0.000000  0.000000
N  161.12676  161.126756   0.000000   0.000000   0.000000 0.000000  0.000000
M -185.02694 -185.026937   0.000000   0.000000   0.000000 0.000000  0.000000
L  110.35247    0.000000   0.000000   0.000000   0.000000 0.000000  0.000000
K    0.00000  161.126756   0.000000   0.000000   0.000000 0.000000  0.000000
J  161.12676    0.000000   0.000000   0.000000   0.000000 0.000000  5.249698
D    0.00000    0.000000   0.000000 180.288684 102.282984 4.259948 -8.534637
C    0.00000    0.000000 180.288684   0.000000  43.301484 1.803448 -3.613137
A    0.00000    0.000000 102.282984  43.301484   0.000000 1.023148  0.000000
B    0.00000    0.000000   4.259948   1.803448   1.023148 0.000000  0.000000
F    0.00000    5.249698  -8.534637  -3.613137   0.000000 0.000000  0.000000
attr(,"call")
nb2mat(neighbours = w, style = "B")
sum(pmw==0)
[1] 132

Pruebas de significancia estadística para el I’ Moran

Ahora podemos probar la significancia. Primero analíticamente, usando lógica basada en regresión lineal y suposiciones.

moran.test(p$hog_NBI_p, ww, randomisation=FALSE)

    Moran I test under normality

data:  p$hog_NBI_p  
weights: ww    

Moran I statistic standard deviate = 1.7582, p-value = 0.03935
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
       0.16093780       -0.07692308        0.01830159

En lugar del enfoque anterior, debe usar la simulación de Monte Carlo. Ese es el método preferido (de hecho, el único método). Funciona de manera que los valores se asignan aleatoriamente a los polígonos, y se calcula el I de Moran. Esto se repite varias veces para establecer una distribución de valores esperados. El valor observado de I de Moran se compara con la distribución simulada para ver qué tan probable es que los valores observados puedan considerarse aleatorios.

moran.mc(p$hog_NBI_p, ww, nsim=99)

    Monte-Carlo simulation of Moran I

data:  p$hog_NBI_p 
weights: ww  
number of simulations + 1: 100 

statistic = 0.16094, observed rank = 97, p-value = 0.03
alternative hypothesis: greater

Diagrama de dispersión de I’ Moran

Se puede hacer un “diagrama de dispersión de Moran” para visualizar la autocorrelación espacial. Primero obtenemos los valores vecinos para cada valor.

n <- length(p)
ms <- cbind(id=rep(1:n, each=n), y=rep(y, each=n), value=as.vector(wm * y))

Se remueven los valor iguales a 0.

ms <- ms[ms[,3] > 0, ]

Se calculan los valores de los promedios de los vecinos.

ams <- aggregate(ms[,2:3], list(ms[,1]), FUN=mean)
ams <- ams[,-1]
colnames(ams) <- c('y', 'y espacialmente retrasado')
head(ams)
      y y espacialmente retrasado
1 17.65                  10.68000
2  7.69                  11.31000
3 18.75                  14.96500
4  4.76                   4.00000
5  0.00                  12.01000
6  4.00                  10.03333

Para finalmente graficar

plot(ams)
reg <- lm(ams[,2] ~ ams[,1])
abline(reg, lwd=2)
abline(h=mean(ams[,2]), lt=2)
abline(v=ybar, lt=2)

Tenga en cuenta que la pendiente de la línea de regresión:

coefficients(reg)[2]
 ams[, 1] 
0.1905849

es casi el mismo valor que el I. de Moran

Aquí hay un enfoque más directo para lograr lo mismo (pero es de esperar que lo anterior aclare cómo se calcula esto en realidad). Tenga en cuenta la estandarización de filas de la matriz de ponderaciones:

rwm <- mat2listw(wm, style='W')
# Revisa si las filas suman hasta 1
mat <- listw2mat(rwm)
apply(mat, 1, sum)[1:15]
   H    I    G    O    N    M    L    K    J    D    C    A    B    F <NA> 
   1    1    1    1    1    1    1    1    1    1    1    1    1    1   NA

Se vuelve a graficar

moran.plot(y, rwm)

Analice el gráfico correspondiente

  1. ¿Cuál es su opinión con respecto a este método?

  2. ¿Tiene importancia el análisis espacial de fenómenos mediante el I’ Moran?

Referencias

Barrantes-Sotela, O., & Solano-Mayorga, M. A. (2020). Comportamiento espacial de la pandemia COVID-19 en Costa Rica durante los meses de marzo y abril de 2020 mediante un análisis de autocorrelación espacial. 1-17.